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ABSTRACT 

We present an analytical calculation of the extreme value statistics for dark matter 
halos - that is, the probability distribution of the most massive halo within some 
region of the universe of specified shape and size. Our calculation makes use of the 
counts-in-cells formalism for the correlation functions, and the halo bias derived from 
the Sheth-Tormen mass function. 

We demonstrate the power of the method on spherical regions, comparing the 
results to measurements in a large cosmological dark matter simulation and achieving 
good agreement. Particularly good fits are obtained for the most likely value of the 
maximum mass and for the high-mass tail of the distribution, relevant in constraining 
cosmologies by observations of most massive clusters. 
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1 INTRODUCTION 

Extreme value (or Gumbel) statistics are concerned with 
the extrema of samples drawn from random distributions. 
If a large number of samples are drawn from some distribu- 
tion, the Central Limit Theorem states that their respective 
means will follow a distribution which tends, in the limit of 
large sample size, to a member of the family of normal dis- 
tributions. Analogously, the maximum (or minimum) values 
u of the samples will have a distribution whose large sample 
size limit - where such a stable limit exists Q - is a member 
of the family of Gen eralised Extrem e Value (GEV) distribu- 
tions as detailed bv lGumbell l|l958f ): 



■ lnPGEv(i/) = (1 + 7J/) 



-1/7 



y={u~a)/l3. (1) 



The shape parameter 7 is sensitive to the underlying distri- 
bution from which the maxima are drawn, while a and /3 
are position and scale parameters. 

Despite their wide use in other fields, extreme 
value statistics have historically seen very little applica- 
tion in astrophysics; some exceptions are the work of 
iBh avsar fc Barrow, (^^985^ on the brightest galaxies in clus- 
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^ Although it can be shown that where a stable limiting distri- 
bution exists it will take the form JTJ, certain pathological dis- 
tributions give no such limit. For our purposes it is sufficient to 
note that the limit indeed exists for distributions which are of ex- 
ponential type, meaning the cumulative distribution function F 
obeys limx^ao d/dx{{l — Fix))/ F' (x)} = 0, and that this class 
includes all physical distributions relevent to our applications. 



ters and the study of lColesI I 19881) on the hottest hot spots of 
the cosmic microwave background temperature fluctuations. 

The past year or so, however, has shown a resurgence 
of interest in the application of extreme value statistics 
to cosmology and questions of extreme structures, as re- 
vealed either in the clusteri ng of galaxies (|Antal et al.ll2009l : 
lYamil a Yarvura et al."2010'). th e prevalence of mass ive clus- 
ters polz & Pcrlmuttor 20 1^; ICa.v6n et all I2OI0I) or t he 
temperature extrema of the CMB i Mikelsons et al.ll2009l '). 

In this paper, we are interested in the dark matter halos 
of massive galaxy clusters. The number density of extremely 
massive clusters is indeed a sensitive probe of the effects of 
the underlyi ng cosmological m odel and laws of physics on 
large scales (|Mantz et al.ll201ol '). These include for instanc e 
the equation of state of dark ene rgy (iMantz et al.l |2008| ) . 
the possibility of modified gravity (iRapetti et al. 2010l ) the 
physical properties of neutrinos (|Mantz et al. 2010l) and 
prim ordial non-Gaussian density fluctuations ( Cavon et al.l 
I2OI0I ). Although the majority of the previously mentioned 
studies focus on the growth of massive clusters, simply 
knowing the mass of the most massive cluster in a sur- 
vey can already provide strong constraints on cosmology 
||Ho1z fc Perlmutterllioioh . 

In the present work, we outline an analytical derivation 
of the extremal halo mass distribution in standard cosmolo- 
gies with Gaussian initial conditions. Rather than taking a 
phenomenological approach, we aim to predict the distri- 
bution of the most massive halo in a region for any speci- 
fled combination of power spectrum, cosmological parame- 
ters and region size and shape. The paper is organized as 
follows. In § [5] we outline the basics of our method for ob- 
taining an analytical expression of the Gumbel distribution 
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of most massive clusters masses and make an explicit link 
with eq. ((T| . In § (3] the theoretical predictions are checked 
against measurements in a very large A^-body cosmological 
simulation. Finally, §|4] follows with a short summary of the 
main results and conclusions. 



2 THEORY 

Consider a large patch of the universe, which can be thought 
of as representing the space covered by a volume-limited 
sample of clusters, and denote by mmax the mass of the most 
massive dark matter halo in that patch. We wish to study 
analytically the Gumbel statistics, that is the probability 
distribution function PG(wmax)dminax of the values taken 
by rrimax if we sample a large number of such patches. Ob- 
viously, this distribution will depend on the size and shape 
of the patch, as well as its redshift. 



2.1 General expression of the Gumbel statistics 

Let us define the cumulative Gumbel distribution by 

PG(m) = Prob.(mniax ^ m) = / PG(mnmx)dmmax. (2) 



Such a probability is also the probability Po{m) that the 
patch is empty of h alos of mass above the threshold m 
l|Colombi et alJbOlOl '). hence 



, . dPo 



(3) 



Note that this assumes that there are no significant edge 
effects, i.e. that the boundaries of the catalog do not cross 
too many clusters. This effect is negligible if the patch is 
large compared to the halo size (and sufficiently compact). 

If halos are unclustered then the void probability follows 
simply from Poisson statistics. 



Po{m) = exp(— nV), 



(4) 



where V is the volume of the patch and n = n(> m) 
the mean density of halos above mass m, with the 
appropriate spatial average with redshift made implicit 
n(> m) = (n[> m, z{x)])a:ev- 

We expect the Poisson limit to be reached for patchs of 
size above a few hundred Mpc, where the matter distribu- 
tion becomes homogeneous. Below this patch size, however, 
halos are significantly clustered. In that case, the calculation 
of the void probability can be performed using a standard 
count-in-cell formalism if the connected A'^-point correla- 
tions functions, ^%(xi, . . . , xn), of halos above the threshold 
are k nown (e.g.. ISzapudi fc SzalavllToQ^ : iBalian fc Schaeffeij 
Il989l ). The superscript h in the previous expression indicates 
halo correlation functions, while the naked ^jv refer to cor- 
relations of the underlying matter density field. Since devi- 
ations from Poisson behavior occur only for moderate patch 
sizes, the complex lightcone effects on the correlations in- 
duced by the evolution of clustering with redshift inside the 
patch (e.g. iMatsubara et al., 199 7') can safely be neglected. 

In particular, one can define the averaged correlations 
over a patch of volume V: 



1 

yW 



■d^XN^%{xi, 



,Xn) 



(5) 



and the typical number of halos above the threshold m in 
excess to the average in overdense patches as: 

TVc = n1/g. (6) 

In the high-m limit, the void probability can be written 

Po(m) =exp[-nV<7(iVc)], (7) 

where the function (t(j/) reads 

1 



l+-f^le 



= y, 



(8) 



(jBernardeau fc Scha effer"l999^1. Note that, as pointed out by 
these authors, this expression for (j{y) follows from a spe- 
cific hierarchical behavior of higher-order correlation func- 
tions of very massive halos at large separations, ^ ~ 

We now proceed to specify the cumulative halo number 
density and the average two-point correlation function of 
halos in order to fully determine the Gumbel statistics. 



2.2 Halo number density 

The number density n(m, z) of halos at a given mass m and 
redshift z, a. k.a. the halo mass functi on, we adopt is the one 
calculated bv lSheth fc TormenI lll999D. It is based on a mod- 
ification of the original model of Press fc Schechted (|l974l ). 
which links the statistics of the initial matter density field 
to the distribution of virialized dark matter halos through a 
spherical top hat description of their gravitational collapse. 
As a result, this mass function can be expressed as a uni- 
versal function of = {dc/u{m, z))'^ , where a{m,z)^ is the 
variance of the initial density field smoothed over spheres of 
radius R{m) containing an average average mass m linearly 
extrapolated to z, and Sc is the critical overdensity threshold 
needed to turn an initial spherical top hat density pertur- 
bation into a collapsed halo at redshift z. More specifically, 
the number density of mass-m halos is given by: 

m^!^^J?i^=A(l + (a.)-)(|^)''%— (9) 



d log u 



with pm = fimP the averaged matter density of the Universe. 
The shape of this mass function is parameterised by a and 
p, and A is simply a normalisation factor. 



2.3 Halo correlation functions 

At sufficiently large separations, the two-point correlation of 
halos of mass m can be related to that of the matter density 
through the bias function 



^2{xi,X2,z) = b{m,z)^^2{xi,X2,z), 



(10) 



where ^2 is the linear dark matter density autocorrela- 
tion at redshift of interest. The function b{m, z) can be 
estimated analytically using the Press-Schechter formalism 
(Mo fc White 1996). Here, to remai n consistent with eg. Q , 
we use the expression for the bias o f ISheth fc TormenI (|l999l ) 



6(m, z) 



1 + 



(1 -f avY ■ 



(11) 



This result is valid in the regime where the separation 
x = \x2 — x\\ is large enough compared to the mass scale 
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Rim). This has been t ested successfully against A^-body sim- 
ulatio ns bv lMo et al.l (fTogQ ) (see however e.g. [Tinker et al] 
|2010| . for possible improvements on eg. Ilip . 

We obtain the bias of halos exceeding mass threshold 
m by calculating the weighted average 



6(> m, z) — 



b(m' , z)n{m' , z)dm' 
n{m', z)dm' ' 



(12) 



and hence the averaged two-point correlation function for 
halos above the threshold, 



C2(> m,z) = b(> m,zf^2(z). 



(13) 



Recall that this equation should be valid in the regime where 
the patch size is large compared to R(m), 



(14) 



but small enough that light cone effects on the clustering 
inside it are negligible. 



Sheth-Tormen, best fit 
Sheth-Tormen, original 
Simulation 




M / Ms,„ 



Figure 1. The upper panel shows the mass function of halos 
in the simulation (points), compared to the Sheth-Tormen mass 
function with (p, a) equal to (0.3, 0.707) and (-0.19, 0.777) (solid 
blue and dashed green lines respectively). The lower panel shows 
the residuals of the two theoretical curves compared to the data; 
i.e. (theory — data)/data. 



2.4 Generalised Extreme Value parameterisation 

The method outlined above allows us to compute the com- 
plete distribution function of of the most massive clusters. 
However, due to its complexity and the necessity of com- 
puting some of the integrals numerically, it does not provide 
us with a neat analytic parameterisation of the distribu- 
tion. Therefore, in order to achieve such a parameterisation, 
we turn to the general theory of extremes, eq. ([l|, using 
u — logj^Q m as the random variable. In order to calculate 
the parameters 7, a and 13, we perform Taylor expansions 
of the analytic Pgev, and Po as computed by our method 
in the Poisson regime (eq. |4|. This Taylor expansion is per- 
formed about the peaks of the two distributions dPo/du and 
dPoBv/du. Equating the first two terms in these expansions 
give us expressions for the three parameters: 



7 = n{> mo)V — 1, /? = 



(l-f7)'^+"' 
n{mo)Vmo In 10 



a = logiomo [(1 + 7) 



(15) 



where mo is the mass at which the distribution dPo/dz = 
In(lO) m pG{m) peaks - hence close to the most likely value 
of m - and is given implicitly by 



A 



PmV 

mo 



-aUQ/2 



(1 + (auo)-") = 



2 2uo l + {aiyo)-P v'i}' 



(16) 



where fo, v'^ and v'^ are v and its derivatives with respect 
to m evaluated at m = mo. 

These equations, then, allow us to neatly summarise 
the information contained in the extreme value distribution 
with the single parameter 7 which describes its shape. This 
statistic has the potential to be used as a tool to compare 
model s with data or with each other, as iMikelsons et al.l 
120091 ') proposed for the CMB. 



3 NUMERICAL EXPERIMENT 

To test our halo mass Gumbel distribution we compare the 
analyti cal result to measur ements on the Horizon 4H Simu- 
lation iTevssier et al.|[2009l ). a large cosmological dark mat- 
ter simulation performed using the RAMSES A'^-body code 
(lTevssieijl200i i. The simulation followed the evolution of a 
cubic piece of the universe 2/i~^Gpc on a side containing 
4096'' particles, i.e. with a particle mass 
Initial conditions w ere based on the WMAP 3-year results 
(|Spergel et al.ll2007l ). with the Hubble constant, density and 
characteristic parameters of the power spectrum given by 
(h, Q.K, n^, fib, erg, n,) = (0.73, 0.76, 0.24, 0.042, 0.77, 
0.958). Halos in the simulation were identified at present 
time, jz = 0, using a 'Friends-of-Friend s' algorithm (e.g. 
IZeldovich et al ] 1 19821 : iDavis et al.l Il985l l with a standard 
linking length parameter value given by 0.2 times the mean 
interparticle distance. 



3.1 Fit of the mass function 

Any discrepancies between our derived Gumbel distribu- 
tion and the true distribution can be thought of as arising 
from one of two sources: either inaccuracies in our chosen 
mass function, or inaccuracies due to the various assump- 
tions made in proceeding from the mass function to pQ. 
In order to quantify the respective contributions of each 
of these two sources, we repeat our calculations with two 
sets of parameters for t he mass function : (i) o nce taking 
the parameters used in ISheth fc TormenI () 19991 ). (p, a) = 
(0.3, 0.707), and, since the Sheth & Tormen form is its stan- 
dard parame trisation is known t o perform on ly approxi- 
mately (e.g. .WarTen et al.l [ioO^ : I Jenkins et al.l ,20011 (ii) 
once with a best fit for a and p to the simulation's mass 
function, leading to (p, a) = (—0.19,0.777). For this latter, 
we also weight bins by their mass, since it is the high-mass 
end of the distribution which is of interest to us. Fig.[T]shows 
both these mass functions along with that measured in the 
simulation. 
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mmax/Msun rrimax/Msun 

(a) L = 100/h Mpc; 7 = -0.11 (b) L = 50/h Mpc; 7 = -0.14 



• • ■ Best fit with Poisson 
Best fit witin bias 




mmax/Msun rrimax/Msun 

(c) L = 20/h Mpc; 7 = -0.11 (d) L = 8/h Mpc; 7 = -0.12 



Figure 2. Distribution of largest cluster masses rrimax. The y axis is probability density per log mass. Points are measurements from 
the simulation with Poisson error bars for each mass bin. The solid blue line is the theoretical result using the best-fit mass function 
(p = —0.19, a = 0.777) and full halo clustering. Green dashes have instead the original Sheth-Tormen parameters (p = 0.3, a = 0.707). Not 
surprisingly, they do not agree as well with the measurement as the solid blue line; red dots have the best-fit values but assume halos are 
Poisson distributed. The orange dot-dashed line is the Generalised Extreme Value distribution, with parameters calculated as explained 
in section [2.41 and assuming Poisson statistics. All calculations use spherical patches, with radius L = 100, 50, 20, 8h~^ Mpc respectively. 



3.2 Results 

Fig. [2] shows the distribution pG(?7imax) calculated as above, 
both for Poisson statistics (equation |4j red dots) and incor- 
porating full clustering (equation [T] blue solid lines) using 
our best-fit mass function. The full clustering calculation is 
also shown for the original Sheth-Tormen parameters (green 
dashes). Points show measurements from the Horizon 411 
simulation for comparison. 

Fitting a Gumbel distribution eq. ^ to the data pre- 
sented in the four panels of this figure yields a single value of 
7 around -0.21 ioioi with reduced ^ 1-1, whereas the an- 
alytic prediction presented in section gives —0.14 ^ 7 ^ 
—0.1. This lack of agreement has its root in the fact that 7 
is very sensitive to the higher order (skewness and kurtosis) 
moments of the data distribution around its peak, and these 
are poorly captured by assuming Poisson statistics, even on 
WOh~^ Mpc scales. 

However, for a patch (in this case a sphere of radius L) 
of size L — 100/i~^ Mpc, Fig. 2(a) shows that pG('Timax) 



theoretical results as it would have seemed from the value 
of 7 alone: we are closing in on the so-called "scale of homo- 
geneity" above which the matter distribution is essentially 
unclustered. Reducing the patch size to 50/i"^ Mpc (fig. 



2(b) I causes the full clustering and Poisson curves to diverge. 



As expected, only the calculations incorporating clustering 
remain a good match to the simulation on these smaller 
scales. 



Decreasing L further (fig. 2(c) 2(d) I causes even the 



measured in the data is not as badly described by Poisson 



calculations including clustering to diverge from the data 
as our approximations - in particular the expression for the 
function a{y) in section I^TT] and the condition (|14p - fail out- 
side the large-!/ limit. For instance, we find _R(10^''Mq) — 3 
Mpc//i, and R{W^'^Mq) = 6.4 Mpc/h which is a signifi- 
cant fraction of the respective patch sizes of 8 Mpc//i and 
20 Mpc//i. Despite these limitations, the description of the 
data is still significantly better than that of the Poisson cal- 
culation and remains impressive at the high-mass end. This 
is excellent news as it is this high-mass tail of the distri- 
bution which is of interest for assessing the significance of 
rare events such as surprisingly massive clusters observed in 
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X-ray or redshift surveys. Indeed the lower-mass tail of the 
distribution for which the prediction fails most significantly 
lies at masses below ~ IO'^^Mq, which corresponds to halos 
containing one to a few galaxies rather than tens or hun- 
dreds, and therefore are of limited interest in the search for 
the most massive cluster. 

Moreover, note that the position of the peak of the prob- 
ability distribution function - that is, the most likely value of 
log 

10 JTT-max - is fairly accurately predicted by the theory even 
when the shape of the curve begins to diverge from that of 
the data. Fig. [3]shows this most likely value, log^Q rhmax, as 
a function of L for both theoretical estimates and the sim- 
ulation data (central line and points in the figure). We also 
show in this figure the 95% confidence region on log^g »7imax 
(upper and lower lines and bars). This too is well fit by the 
theory, particularly for the upper limit, and we emphasize 
that this is a crucial test of the theory's ability to give sig- 
nificance values for observations of specific overly massive 
clusters. 

In addition to the four values of L shown in Fig. [2l Fig. 
[3] has a final simulation point at L = 500//i Mpc. Here the 
95% confidence region is poorly fit, because the patch is too 
large compared to the simulation box to get good statistics. 
However, the measured value of log^Q mmax is still in good 
agreement with the theory. 



3.2.1 Senstivity to the mass functton 

Worthy of note is the close similiarity of the curves for the 
two sets of parameters in the Sheth-Tormen mass function 
(solid blue and dashed green lines in Fig. [3)). Although the 
mass functions themselves differ by a relatively large amount 
compared to the best-fit function's agreement with the sim- 
ulation points (Fig. [T] lower panel), this translates to only 
a modest change in the distribution of rrimax. Simil ar calcu- 
lation s per formed with the ma ss functions of .Tinker et al.l 
(|2008l ) and lJenkins et all (|200ll ) lead us to conclude that the 
extreme value distribution is fairly robust to the choice of 
any reasonable analytic fit to the mass function. 



3.2.2 Redshift variation 

While the above calculations use patch sizes L small enough 
that the redshift evolution within the patch is negligible, it 
is also interesting to calculate mmax for a larger region with 
significant Az. As noted in section 12.11 redshift variation 
can be taken into account by a weighted spatial average 
provided the Poisson approximation holds, which is the case 
for such large patches. In particular, averaging the number 
density of halos over the range 2 = to oo gives a value for 
the expected largest mass cluster in the entire observable 
universe. 

We performed this calculation assuming Poisson statis- 
tics, and found mmax = 4.6 ijig x 10^^ Mq at the 1 — a con- 
fidence level. We note that this i s in fair agreement with th e 
similar calculation performed bv lHolz Sz Pe rlmutterl (|201(J ). 
who obtained mmax = 3.8 ±0.5 x 10^^ using a WMAP 7 cos- 
mology and the mass function of [Tinker etall (|2008l ). 




10 100 1000 

L h/Mpc 



Figure 3. The most likely value of logjg mmax (middle line and 
points) and the 95% confidence limits (upper and lower line and 
bars). Dashed red and solid blue lines are analytic results for 
the Poisson limit and full clustering calculations respectively, and 
points are simulation values for L corresponding to the four panels 
of fig. [2] plus L = 500//1 Mpc. 



4 CONCLUSION AND DISCUSSION 

We have presented an analytic prediction of the the prob- 
ability distribution of mmax, the most massive dark mat- 
ter halo/galaxy cluster in a specified region of the universe, 
making use of the counts-in-cells formalism. Our calcula- 
tion, valid for Gaussian initial conditions, is numerically 
consistent with that proposed bv lHolz fc Perlmuttej (12010 ) 
when performed assuming such massive halos are Poisson 
distributed spatially. However, the work presented in this 
paper improves on the calculation performed by these au- 
thors in two aspects: (i) our results are given in a fully ex- 
plicit analytic form and (ii) they include the contribution of 
clustering of halos. 

We also compared our analytic predictions to measure- 
ments from a large (2000/i~^ Mpc on a side) and well re- 
solved (particle mass 7.7h~^ x 10^ M©) cosmological A^-body 
simulation at zero redshift. We achieve remarkable agree- 
ment with the simulation in the area of parameter space in 
which our formalism is expected to be valid, namely patch 
radii above a few tens of h^^ Mpc. More surprisingly, even 
outside this range of scales the high-mass tail of the distribu- 
tion is well fit by our "fully clustered" theoretical estimate, 
as is the most likely value of logjQ mmax. 

This unexpected success over a wide range of scales war- 
rants the application of the formalism to quantify the statis- 
tical significance of individual clusters observed in surveys. 
By applying our method to a patch of shape, size and red- 
shift equivalent to a real survey we can obtain the Gumbel 
distribution and hence a likelihood for the observed value 
of mmax- Moreover, we are quite confident that our method 
can be extended to non-standard cosmologies such as those 
including initial non-Gaussianities. It could therefore pro- 
vide a measure of the evidence for such cosmologies from 
existing surveys of the most massive clusters as advocated in 
ICavon et al](|2010l ). We plan to tackle this exciting prospect 
in the very near future. 

In addition to the full likelihood curve of mmax we are 
able, via the analytic Generalised Extreme Value formal- 
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ism, to produce a summary of the distribution in the form 
of the three parameters a, b and 7. The latter in particular 
has been propos ed previously as a stat istic for use in model 
comparison (e.g. iMikelsons eralll2009l l. Although the work 
described in this paper uses a single cosmology and power 
spectru m and produce s a rou ghly constant value of 7, the re- 
sults of lColombi et al.] (|2010i ) suggest that 7 should be quite 
sensitive to the shape of the power spectrum. If the effect on 
7 of the clustering can be fully quantified, we therefore ex- 
pect to be able to use it as a statistic for direct comparison 
of models with observation. Likewise a, which is closely re- 
lated to the peak mo of the distribution, may prove a useful 
statistic since we have demonstrated that it is well predicted 
by our theory. 
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